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PDE  AND  TENSOR  BASED  APPROACH  FOR  PARIETAL  DYNAMIC 
TRACKING  IN  ULTRASOUND  SEQUENCES 
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Abstract-  A  method  for  pericardial  and  endocardial  edge 
extraction  of  the  left  ventricle  (LV)  in  2-D  short  axis 
echocardiographic  sequence  is  presented.  The  proposed 
algorithm  consists  in  two  steps:  nonlinear  anisotropic  diffusion 
filtering  and  directional  geodesic  active  contours.  These  steps 
are  based  on  both  PDE  and  tensors.  Tensor  formalism  gives 
access  to  structure  orientation  of  the  acoustic  density  field.  This 
provides  a  more  robust  algorithm  by  using  the  whole  image 
information.  Sample  experimental  results  are  shown  for 
myocardium  edge  detection  in  2-D  short-axis  ultrasound 
sequences. 

Keywords  -  echocardiography,  nonlinear  anisotropic  diffusion, 
directional  geodesic  active  contours,  PDE,  tensor. 

I.  INTRODUCTION 

Color  kinesis  is  an  echocardiographic  technique  based  on 
acoustic  densitometry  to  assess  parietal  movements.  It  has 
been  designed  to  make  the  detection  of  contraction 
abnormality  easier  and  is  provided  as  an  online  tool  on 
ultrasound  imaging  system  ( e.g .  Hewlett-Packard  Sonos 
5500).  For  each  frame,  every  pixel  within  a  ROI  outside  the 
LV  cavity  is  classified  either  as  blood  or  endocardium  with 
respect  to  its  integrated  backscatter  level.  Diastolic 
endocardial  motion  is  visualized.  Class  transitions  are 
encoded  using  a  specific  color  hue  for  each  consecutive 
frame,  each  color  corresponding  to  the  wall  excursion  over  a 
33-ms  period  of  time.  Concentric  color  bands  are 
superimposed  on  the  end-diastolic  frame  to  depict  the  entire 
LV  endocardial  excursion  throughout  diastole.  This  technique 
requires  a  high  image  quality  in  order  to  detect  closed  edges. 
To  overcome  this  limitation,  deformable  models  based  on 
PDE  and  tensors  are  used  for  the  tracking  of  myocardium 
walls  from  an  ultrasound  sequence. 

The  present  paper  describes  how  tensor  formalism  can  be 
coupled  to  PDE  in  a  segmentation  problem  in  order  to  better 
take  into  account  gray  level  information  than  with  a  gradient. 

In  section  II,  we  explain  the  relevance  of  tensor  use,  its 
mathematical  framework  and  the  adaptation  of  algorithms  to 
ultrasound  context.  Finally  section  III  focuses  on  results  and 
discussion. 

II.  METHODOLOGY 

Our  method  can  be  divided  into  two  parts.  First,  there  is  a 
preprocessing  step  based  on  nonlinear  anisotropic  diffusion. 
It  aims  at  reducing  speckle  in  the  original  image  without 
altering  its  significant  features.  Secondly,  in  order  to  extract 
myocardium  boundaries,  a  method  of  directional  active 
contours  is  used.  These  processings  are  based  on  both  tensor 
and  PDE. 


A.  The  Structure  Tensor 

In  case  of  low  contrast  or  textured  images,  the  use  of 
simple  edge  detector  like  the  Gaussian  smoothed  gradient 
may  not  be  efficient  enough.  A  more  sophisticated  descriptor 
is  thus  required  by  means  of  structure  tensor. 

ma=v{Ka*i),  (i) 

where  I  is  the  image  gray  level,  K<j  a  Gaussian  kernel  with 
standard  deviation  o  and  *  the  convolution  product.  In  this 
case,  o  is  usually  called  local  space. 

The  previous  regularized  gradient  gives  a  local  direction, 
but  a  more  relevant  information  is  the  orientation.  A  tensor 
product  is  defined  by: 

=  V/„'(V/J  (2) 

This  matrix  owns  interesting  intrinsic  properties.  Indeed, 
its  eigenvectors  (v!  v2)  define  an  orthonormal  basis,  and  the 
corresponding  eigenvalues  (|i|  ,  |i2)  give  a  measure  of 
contrast  at  the  local  scale  o: 

V.//V4  v2±V4 

,  ,2  (3) 

A=iv/,r  ^2=° 

Finally,  the  convolution  of  J  0  (VY^. )  with  a  Gaussian 

kernel  Kp  gives  a  local  descriptor  of  gray  level  orientation 
called  structure  tensor: 

Je(V4)  =  ^p*J„(v/J  (4) 

The  standard  deviation  p  of  this  kernel  is  called  the 
integration  scale.  J  (V/CT)  reflects  the  image  coherence.  A 

geometric  interpretation  can  be  inferred  from  spectral 

decomposition  [11],  Indeed,  the  eigenvector  associated  to  the 
highest  eigenvalue  A.)  corresponds  to  the  average  direction  of 
the  maximal  gray  level  variation.  Furthermore,  the  other 
eigenvector  gives  the  dominant  local  direction,  called 

coherence  direction.  The  eigenvalues  (A,  >  A2  >  0)  describe 
the  local  configurations  of  gray  levels: 

>  constant  areas:  Ai  =  A2  =  0, 

>  edges:  Ai  » A2  =  0, 

>  comers:  Ai  »  A2  »  0. 


Therefore,  Ai  -  A2  acts  as  a  local  measure  of  anisotropy. 
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B.  Non-Linear  Anisotropic  Diffusion 

Diffusion  in  image  processing  is  equivalent  to  a  physical 
process,  which  aims  at  equilibrating  differences  in 
concentration  (Tick’s  law)  under  the  mass  condition 
conservation.  A  mathematical  formulation  leads  to  the 
diffusion  equation: 


d,C  =  div{DVC)  (5) 


where  C  is  a  concentration  and  D  a  diffusion  tensor.  This 
equation  is  the  familiar  heat  equation  in  the  context  of  heat 
transfer. 

In  image  processing,  we  identify  the  concentration  with 
the  image  gray  level  at  a  given  location.  This  equation  is  the 
starting  point  of  filtering  diffusion  methods  [15].  From  all  the 
existing  methods,  three  different  classes  can  be  extracted: 

>  linear  isotropic  diffusion  filter  with  constant  scalar 
diffusivity  [7,  17]; 

>  nonlinear  isotropic  diffusion  filter  with  scalar  diffusivity 
adapted  to  the  local  image  structure  [4,10]; 

>  nonlinear  anisotropic  diffusion  filter  with  diffusion 
tensor  adapted  to  the  local  image  structure  [5,16], 

The  nonlinear  isotropic  diffusion  is  not  suited  to 
ultrasound  image  because  of  noisy  edges.  On  the  other  hand 
the  nonlinear  anisotropic  diffusion  is  able  to  eliminate  noise 
while  enhancing  edges. 

Finally,  the  model  of  nonlinear  anisotropic  diffusion  with 
temporal  regularization  proposed  by  Cottet  and  El-Ayyadi  [5] 
is  given  by  a  coupled  PDE-ODE  system: 

=  div[p  VI ), 

vt  (6) 

~~  =  “  (fs  (v/)  -  d). 

dt  T 


X  is  a  diffusion  parameter  and  Fs  (V/)  a  nonlinear  operator 
of  orthogonal  projection  defined  as: 


F.(v/)=ur1_M\d+jvj 


if  |V/|  >  s 
if  not 


with  PVL1  = 


|v/|2  L 


Iy  7tMvi=(ix  ii 

rjy  n  1  v  y 


(7) 


where  s  is  a  contrast  parameter  which  is  a  threshold  on  image 
gradient  selecting  significant  edges. 

Operator  PVL1  can  be  interpreted  as  the  +7t/2  rotation  of 
n  means  that  image  is  diffused  along  edges 

perpendicularly  to  the  orientation  defined  by  the  structure 
tensor. 


C.  Directional  Geodesics  Actives  Contours 


The  model  of  geodesic  active  contour  was  introduced  by 
Caselles  [3]  as  a  geometric  alternative  for  snakes,  which  are 
the  parametric  active  contour  of  Kass  et  al  [8]. 

Snakes  are  based  on  an  energy  minimization  along  a 
parametric  curve.  This  energy  is  composed  of  an  internal 
term  and  an  external  one.  Two  major  problems  may  be 
encountered.  First,  the  energy  non-convexity  does  not 
guarantee  the  uniqueness  of  the  solution.  Secondly,  topology 
changes  and  singularities  are  forbidden  by  the  explicit 
representation. 

The  geodesic  active  contours  [3]  are  not  only  based  on  an 
energy  minimization,  but  are  also  parameter-free  thanks  to  an 
intrinsic  geometric  model.  Evolution  equation  uses  the  level 
set  representation  of  Osher  and  Sethian  [9]: 


f=Hv- 


V  v  "  J 


f =s(Kl)N 


Vv“' 


+lSl(K|)|V„| 
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In  this  relation  u  is  the  implicit  contour  representation.  In 
the  level  set  formulation,  inward  unit  normal  vector  of  the 

contour  is  intrinsically  given  by  Vw/|Vw|  and  curvature  by 

V.(Vm/|Vm|).  Concerning  boundary  detection,  the  stopping 

function  g  is  a  positive  decreasing  function  of  the  image 
gradient.  The  right  hand  side  of  (8)  can  be  interpreted  as  the 
contribution  of  two  forces.  The  first  one  is  linked  to  the 
contour  geometry,  weighted  by  the  stopping  function  and 
combines  a  local  contribution  (curvature)  and  a  global  speed 
one  (inner  pressure  parameter  v).  The  second  one  is  an 
attraction  term  oriented  towards  edges. 

Like  snakes,  geodesic  active  contours  use  only  minimal 
geometric  information  by  means  of  a  first  order  isotropic 

estimator  |V/a  |  associated  to  the  stopping  function  g:  it 

weights  the  contour  normal  speed  and  is  part  of  the  attraction 
term. 

A  model  of  directional  geodesic  active  contour  based  on 
tensor  formalism  is  proposed  in  [11]  in  order  to  manage  the 
orientation  of  local  structures.  This  model  can  be  interpreted 
as  an  anisotropic  generalization  of  the  geodesic  actives 
contours.  The  first  model  of  geodesic  active  contour  [3]  is 
generalized  by: 


du 

a7 


=  |Vm|V  ■ 


D 


Vm 

H 


-V'Vdet(D)Vn 


with  the  following  boundary  condition: 

D  Vwit  =  0  on  x  (0,  T  ], 


(9) 


(10) 


where  Q  is  the  image  domain.  In  this  PDE,  D  is  an 
interaction  tensor  computed  as  a  function  of  the  image 
structure  tensor. 
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In  order  to  reflect  the  image  local  geometric  structure,  the 
eigenvectors  of  D  have  to  coincide  with  those  of  Jg  (V/CT  ). 

Associated  eigenvalues  are  replaced  by  g(A,rA,2)  and  1 
respectively  to  specialize  contour  behavior  along  the  image 
structure.  The  interaction  tensor  D  is  a  generalization  of  the 
edge  detector  g. 

D.  Application  to  parietal  tracking  in  ultrasound  sequences 

For  myocardium  segmentation  on  short-axis  view, 
different  approaches  are  to  be  used  for  pericardium  and 
endocardium.  Indeed,  there  are  physiologically  important 
differences  between  these  two  objects.  First,  endocardial  wall 
is  an  interface  between  blood  and  muscle  with  rather  high 
contrast,  while  pericardium  presents  a  contrast  inversion  from 
anterior  walls  to  posterior  ones.  Secondly,  while  endocardium 
has  a  complex  geometry  and  goes  through  large  topological 
changes,  pericardium  shape  is  less  variable  throughout  the 
cardiac  cycle. 

•  Pericardium: 

After  a  manual  initialization  on  the  first  frame,  active 
contour  model  is  used  with  a  zero  pressure  parameter. 
Indeed,  a  small  movement  assumption  is  taken  into  account, 
for  this  boundary,  i.e.  contour  should  stay  in  the  attraction 
valley  of  the  stopping  function  between  two  consecutive 
frames. 

•  Endocardium: 

For  each  frame,  a  seed  is  put  in  the  middle  of  the  image 
and  pressure  parameter  v  is  used  in  order  to  push  the  contour 
towards  the  endocardium.  But  a  constant  parameter  is  not 
relevant,  because  residual  pressure  forces  could  push  the 
contour  through  low  edges.  A  new  parameter  depending  on 
the  global  contour  energy  is  thus  introduced: 

<jy  det(D(ir'  (0)(s)))cfc 

v(ir1(o)(,s'))=  sign(Al).— — - i - ,  (H) 

¥s 

u-'(0) 

where  s  is  the  curvilinear  abscissa  on  the  curve  u  ‘(0),  A  is  the 
laplacian  operator  and  sign,  the  sign  operator. 

The  global  energy  criterion  insures  that  the  contour  slows 
down  even  on  weakly  constrained  segments.  Furthermore,  the 
sign  of  the  image  laplacian  is  applied  to  define  a  return  force 
oriented  towards  edges. 

The  relative  stability  of  |v|  is  used  as  a  stopping  criterion. 

In  practice,  the  numerical  scheme  proposed  in  [5]  is  used 
for  the  diffusion  process.  For  active  contours,  the  implicit 
formulation  of  level-set  [9]  is  used  combined  to  a  narrow- 
band  optimization  [1], 

III.  RESULTS  AND  DISCUSSION 

The  previous  algorithms  are  applied  to  IBS  (Integrated 
BackScatter)  ultrasound  sequences  on  short-axis  view  from  a 
Hewlett-Packard  Sonos  5500  imaging  system.  This 


acquisition  mode  is  usually  used  for  quantification  purposes 
of  cyclic  variations  (CVs)  and  provides  smoother  images. 
Images  are  encoded  as  512*480  matrices  on  a  8-bit  depth  and 
sequences  sampled  at  a  25  Hz  rate. 

A.  Nonlinear  anisotropy  diffusion 

Figure  1  gives  the  result  of  the  nonlinear  anisotropic 
diffusion  filtering.  Three  parameters  are  to  be  adjusted. 

Contrast  parameter  is  chosen  to  preserve  significant 
edges:  a  value  of  0.25  provides  a  good  enhancement  of  the 
endocardial  wall.  The  pericardial  wall  exhibits  a  lower 
contrast  not  far  from  the  noise  level  and  is  thus  less 
pronounced. 

The  diffusion  parameter  x  is  directly  related  to  the 
iteration  number  n.  Indeed,  it  has  been  proved  that  the 
algorithm  tends  towards  a  steady  state  within  a  few  x.  A  good 
compromise  between  diffusion  quality  and  computational 
cost  is  obtained  for  X  =  40  and  n  =  3x. 

B.  Directional  geodesic  actives  contours 

Figure  2  gives  an  example  of  endocardium  detection.  The 
same  process  is  performed  for  each  frame  and  the  result  of 
the  endocardium  tracking  is  shown  for  the  systolic  period. 
Three  contours  out  of  seven  have  been  sampled  for  better 
visualization. 

Concerning  pericardium  tracking,  four  consecutive 
pericardium  walls  have  been  drawn  on  the  same  image  in 
order  to  assess  its  movement  during  the  diastolic  period. 

Small  values  of  the  scale  parameters  (o=0.5,  p=0.5)  are 
chosen  for  these  two  tracking  steps  because  of  the  diffusion 
preprocessing.  In  case  of  higher  values,  corresponding 
Gaussian  filtering  would  alter  edge  information.  The  new 
pressure  parameter  V  (1 1)  increases  convergence  rate. 

Whereas  endocardial  wall  tracking  is  fully  automatic, 
pericardial  wall  one  requires  user  interaction  when  the 
contour  is  too  far  from  the  attraction  valley  ( 1  image  out  of  5 
on  average). 

V.  CONCLUSION 

A  parietal  dynamic  tracking  for  ultrasound  sequences  was 
presented,  which  combines  a  preprocessing  step  based  on  a 
nonlinear  anisotropic  diffusion  filtering  and  directional 

geodesic  active  contours.  A  2-D  version  of  the  algorithms 
was  presented  and  a  3-D  (2-D+T)  generalization  is  planned 
for  nonlinear  anisotropic  filtering.  At  present,  this  method 
suffers  from  high  calculation  cost  and  algorithms  for  both 
nonlinear  anisotropic  diffusion  filtering  and  directional 

geodesic  active  contour  could  be  improved  by  using  a  fast 
algorithm  based  on  parallel  implementations  of  AOS 
(Additive  Operator  Splitting).  [6,14] 

ACKNOWLEDGMENT 

The  authors  wish  to  thank  Pr.  Cassagnes,  Pr.  Lusson,  Dr. 
Motreff  and  Dr.  Langlade  from  the  department  of  cardiology 
of  the  Gabriel  Montpied  hospital  (Clermont-Ferrand,  France) 
for  their  support  and  valuable  suggestions. 


4  of  4 


Figure  1  Pre-treatment  filtering 


Figure  2  Detection  and  tracking  of  endocardium 


Figure  3  Detection  and  tracking  of  pericardium 
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